import numpy as np

movimento_oscilador_forcado
calcular_amplitude_stationaria
calcular_periodo_stationario
calcular_energia_mecanica
energia_dissipada_amortecimento


#P11 - 1
#a)numericamente a lei do movimento

#b)amplitude do movimento e o seu período no regime estacionário, ultimos 20% dos dados (regime estacionário)

#c)valores de 𝜔𝑓 entre 0.2 e 2 rad/s. gráfico da amplitude em regime estacionário em função de 𝜔𝑓 . 
# Qual a frequência angular 𝜔𝑓 que corresponde à maior amplitude?

def movimento_oscilador_forcado(m, k, b, F0, wf, x0, v0, t_max=100, n_pts=2000):
    w0 = np.sqrt(k/m)
    beta = b/(2*m)
    wa = np.sqrt(w0**2 - beta**2)

    # Amplitude e fase estacionária
    A = (F0/m) / np.sqrt((wf**2 - w0**2)**2 + (b*wf/m)**2)
    alpha = np.arctan2(b*wf/m, w0**2 - wf**2)

    # Parte transiente
    phi = np.arctan2(-(v0 + x0*beta)/wa, x0)
    Aa = x0 / np.cos(phi)

    # Solução completa
    t = np.linspace(0, t_max, n_pts)
    x_t = (Aa * np.exp(-beta * t) * np.cos(wa * t + phi) +
            A * np.cos(wf * t + alpha))

    return t, x_t


def calcular_amplitude_stationaria(t, x_t):
    """
    Calcula a amplitude do movimento no regime estacionário.
    """
    N = int(0.8 * len(t))
    x_stationary = x_t[N:]
    amplitude = (np.max(x_stationary) - np.min(x_stationary)) / 2
    return amplitude


def calcular_periodo_stationario(t, x_t, wf):
    """
    Calcula o período do movimento no regime estacionário.
    """
    N = int(0.8 * len(t))
    t_stationary = t[N:]
    x_stationary = x_t[N:]

    T_esperado = 2 * np.pi / wf
    imax1 = np.argmax(x_stationary)
    tmax1 = t_stationary[imax1]

    idx_search = np.where(t_stationary > tmax1 + 0.5 * T_esperado)[0]

    if len(idx_search) > 0:
        imax2 = np.argmax(x_stationary[idx_search])
        tmax2 = t_stationary[idx_search][imax2]
        periodo = tmax2 - tmax1
        return periodo
    else:
        return np.nan
    


def calcular_energia_mecanica(m, k, t, x_t):
    """
    Calcula a energia mecânica total do sistema ao longo do tempo.
    
    Parâmetros:
    m  - massa do oscilador (kg)
    k  - constante elástica da mola (N/m)
    t  - vetor de tempo (s)
    x_t - vetor de posições x(t) ao longo do tempo (m)
    
    Retorna:
    Em - vetor da energia mecânica total em cada instante de tempo (J) (array)
    """
    # Derivada numérica para v(t)
    dt = t[1] - t[0]               # Intervalo de tempo entre amostras
    v_t = np.gradient(x_t, dt)    # Derivada de x(t) para obter v(t) ≈ dx/dt

    # Energia cinética (Ec) e energia potencial (Ep)
    Ec = 0.5 * m * v_t**2          # Energia cinética: Ec = ½ m v²
    Ep = 0.5 * k * x_t**2          # Energia potencial elástica: Ep = ½ k x²

    Em = Ec + Ep                   # Energia mecânica total: Em = Ec + Ep
    return Em

t_i = 37.5  # tempo desejado em segundos
index = np.searchsorted(t, t_i)
Em_ti = Em[index]

print(f"Em({t_i:.1f} s) = {Em_ti:.4f} J")

delta_Em = Em[-1] - Em[0]
print(f"Variação total da energia mecânica: {delta_Em:.4f} J")


def energia_dissipada_amortecimento(b, t, x_t):
    dt = t[1] - t[0]
    v_t = np.gradient(x_t, dt)
    dissipacao = b * np.sum(v_t**2) * dt  # ∫ b*v² dt
    return dissipacao